function bayesian_svm_training()
% 使用贝叶斯优化的高效SVM训练脚本
% 相比传统网格搜索，贝叶斯优化能更快找到更优参数

fprintf('===== 贝叶斯优化SVM训练系统 =====\n');
fprintf('优势：\n');
fprintf('1. ✓ 贝叶斯优化 - 比网格搜索更高效\n');
fprintf('2. ✓ 更大的参数搜索空间\n');
fprintf('3. ✓ 自适应参数调整\n');
fprintf('4. ✓ 更好的收敛性能\n');
fprintf('5. ✓ 支持多核并行加速\n\n');

% 设置随机种子，确保结果可重现
rng(42, 'twister');
fprintf('✓ 已设置固定随机种子，确保结果可重现\n');

% 配置参数
config = struct();
config.trainFolder = 'D:\文档资料\OneDrive\桌面\人脸配准(眉毛特征点)\train';
config.trainImgPerPerson = 5;
config.pcaRatio = 0.95;
config.imgSize = 100;

% 贝叶斯优化配置
config.bayesian = struct();
config.bayesian.maxEvaluations = 50;        % 最大评估次数（比网格搜索20次更多）
config.bayesian.maxTime = 300;              % 最大时间限制（秒）
config.bayesian.useParallel = false;        % 是否使用并行计算
config.bayesian.acquisitionFunction = 'expected-improvement-plus';
config.bayesian.kfold = 5;                  % 交叉验证折数

fprintf('训练配置:\n');
fprintf('- 训练文件夹: %s\n', config.trainFolder);
fprintf('- 每人训练图像数: %d\n', config.trainImgPerPerson);
fprintf('- PCA贡献率: %.2f\n', config.pcaRatio);
fprintf('- 贝叶斯优化最大评估次数: %d\n', config.bayesian.maxEvaluations);
fprintf('- 贝叶斯优化时间限制: %d秒\n', config.bayesian.maxTime);
fprintf('\n');

% 检查训练文件夹
if ~exist(config.trainFolder, 'dir')
    error('训练文件夹不存在: %s', config.trainFolder);
end

% 检查并行计算工具箱
if config.bayesian.useParallel
    try
        parpool('local');
        fprintf('✓ 已启用并行计算加速\n');
    catch
        fprintf('⚠ 并行计算工具箱不可用，使用串行计算\n');
        config.bayesian.useParallel = false;
    end
end

% 1. 数据加载
fprintf('【步骤1/6】加载训练数据...\n');
[trainFaces, trainLabels, fileInfo] = loadTrainingDataStable(config);
fprintf('✓ 加载完成: %d个人，共%d张图像\n', length(unique(trainLabels)), length(trainLabels));

% 2. 数据预处理
fprintf('【步骤2/6】数据预处理和标准化...\n');
[trainFaces, meanFace] = preprocessDataStable(trainFaces);
fprintf('✓ 数据标准化完成\n');

% 3. PCA降维
fprintf('【步骤3/6】PCA降维...\n');
[eigenFaces, projectMatrix, weightedTrain] = performStablePCA(trainFaces, meanFace, config.pcaRatio);
fprintf('✓ PCA降维完成: %d -> %d 维\n', size(trainFaces, 1), size(eigenFaces, 2));

% 4. 贝叶斯优化SVM训练
fprintf('【步骤4/6】贝叶斯优化SVM训练...\n');
startTime = tic;
[svmModel, bayesianResults] = trainBayesianSVM(weightedTrain', trainLabels, config.bayesian);
optimizationTime = toc(startTime);
fprintf('✓ 贝叶斯优化完成，耗时: %.2f秒\n', optimizationTime);

% 5. 模型评估
fprintf('【步骤5/6】模型性能评估...\n');
trainAccuracy = testSVMAccuracy(svmModel, weightedTrain', trainLabels);
fprintf('✓ 训练集准确率: %.2f%%\n', trainAccuracy);

% 6. 独立交叉验证
fprintf('【步骤6/6】独立交叉验证...\n');
cvAccuracy = performIndependentCV(weightedTrain', trainLabels, bayesianResults.bestParams);
fprintf('✓ 独立交叉验证准确率: %.2f%%\n', cvAccuracy);

% 保存贝叶斯优化模型
fprintf('\n保存贝叶斯优化模型...\n');
bayesianModel = struct();
bayesianModel.svmModel = svmModel;
bayesianModel.eigenFaces = eigenFaces;
bayesianModel.projectMatrix = projectMatrix;
bayesianModel.meanFace = meanFace;
bayesianModel.trainLabels = trainLabels;
bayesianModel.config = config;
bayesianModel.trainAccuracy = trainAccuracy;
bayesianModel.cvAccuracy = cvAccuracy;
bayesianModel.bayesianResults = bayesianResults;
bayesianModel.optimizationTime = optimizationTime;
bayesianModel.trainingTimestamp = datestr(now);

save('bayesian_svm_model.mat', 'bayesianModel');
fprintf('✓ 模型已保存到: bayesian_svm_model.mat\n');

% 显示优化结果
fprintf('\n===== 贝叶斯优化结果 =====\n');
fprintf('最优参数:\n');
fprintf('- KernelScale: %.4f\n', bayesianResults.bestParams.KernelScale);
fprintf('- BoxConstraint: %.4f\n', bayesianResults.bestParams.BoxConstraint);
fprintf('- 最优交叉验证损失: %.6f\n', bayesianResults.bestObjective);
fprintf('- 总评估次数: %d\n', bayesianResults.numEvaluations);
fprintf('- 优化时间: %.2f秒\n', optimizationTime);

% 性能对比
fprintf('\n===== 性能评估 =====\n');
fprintf('训练集准确率: %.2f%%\n', trainAccuracy);
fprintf('交叉验证准确率: %.2f%%\n', cvAccuracy);

if trainAccuracy >= 99.0 && cvAccuracy >= 95.0
    fprintf('\n🎉 优秀！模型达到高性能要求！\n');
    fprintf('   训练准确率: %.2f%% (≥99%%)\n', trainAccuracy);
    fprintf('   交叉验证准确率: %.2f%% (≥95%%)\n', cvAccuracy);
else
    fprintf('\n📊 模型性能总结\n');
    fprintf('   训练准确率: %.2f%%\n', trainAccuracy);
    fprintf('   交叉验证准确率: %.2f%%\n', cvAccuracy);
    if trainAccuracy < 99.0
        fprintf('   建议: 考虑增加训练数据或调整PCA参数\n');
    end
    if cvAccuracy < 95.0
        fprintf('   建议: 模型可能存在过拟合，考虑增加正则化\n');
    end
end

% 显示贝叶斯优化优势
fprintf('\n===== 贝叶斯优化优势总结 =====\n');
fprintf('1. ✓ 智能参数搜索 - 不需要预定义网格\n');
fprintf('2. ✓ 高效收敛 - 用更少评估找到更好参数\n');
fprintf('3. ✓ 连续参数空间 - 不限于离散值\n');
fprintf('4. ✓ 自适应策略 - 根据历史结果调整搜索方向\n');
fprintf('5. ✓ 时间控制 - 可设置最大时间限制\n');

fprintf('\n===== 贝叶斯SVM训练完成 =====\n');

end

function [trainFaces, trainLabels, fileInfo] = loadTrainingDataStable(config)
% 稳定加载训练数据，确保每次加载顺序一致

% 获取所有图像文件并排序
imgExtensions = {'.jpg', '.png', '.bmp', '.jpeg', '.tif'};
allFiles = [];

for i = 1:length(imgExtensions)
    pattern = ['*' imgExtensions{i}];
    files = dir(fullfile(config.trainFolder, pattern));
    allFiles = [allFiles; files];
end

if isempty(allFiles)
    error('训练文件夹中未找到图像文件');
end

% 按文件名排序，确保加载顺序一致
[~, sortIdx] = sort({allFiles.name});
allFiles = allFiles(sortIdx);

% 解析文件名并按人物ID分组
personData = containers.Map('KeyType', 'int32', 'ValueType', 'any');

for i = 1:length(allFiles)
    filename = allFiles(i).name;
    tokens = regexp(filename, '(\d+)[-_](\d+)', 'tokens');
    
    if ~isempty(tokens)
        personID = str2double(tokens{1}{1});
        imageSeq = str2double(tokens{1}{2});
        
        if isKey(personData, personID)
            data = personData(personID);
            data.sequences(end+1) = imageSeq;
            data.files(end+1) = i;
            personData(personID) = data;
        else
            data = struct();
            data.sequences = imageSeq;
            data.files = i;
            personData(personID) = data;
        end
    end
end

% 按人物ID排序并稳定选择训练图像
personIDs = sort(cell2mat(keys(personData)));
trainFaces = [];
trainLabels = [];
fileInfo = [];

for i = 1:length(personIDs)
    personID = personIDs(i);
    data = personData(personID);
    
    % 按序号排序，确保选择的一致性
    [sortedSeqs, sortIdx] = sort(data.sequences);
    sortedFiles = data.files(sortIdx);
    
    % 选择前N张图像
    numToUse = min(config.trainImgPerPerson, length(sortedFiles));
    selectedFiles = sortedFiles(1:numToUse);
    
    % 加载图像
    for j = 1:length(selectedFiles)
        fileIdx = selectedFiles(j);
        imgPath = fullfile(config.trainFolder, allFiles(fileIdx).name);
        
        img = imread(imgPath);
        if size(img, 3) == 3
            img = rgb2gray(img);
        end
        
        % 标准化图像尺寸
        img = imresize(img, [config.imgSize, config.imgSize]);
        imgVector = double(img(:));
        
        trainFaces = [trainFaces, imgVector];
        trainLabels = [trainLabels; personID];
        fileInfo = [fileInfo; struct('personID', personID, 'filename', allFiles(fileIdx).name)];
    end
end

end

function [trainFaces, meanFace] = preprocessDataStable(trainFaces)
% 稳定的数据预处理

% 计算均值并中心化
meanFace = mean(trainFaces, 2);
trainFaces = trainFaces - repmat(meanFace, 1, size(trainFaces, 2));

% L2标准化
for i = 1:size(trainFaces, 2)
    norm_val = norm(trainFaces(:, i));
    if norm_val > 0
        trainFaces(:, i) = trainFaces(:, i) / norm_val;
    end
end

end

function [eigenFaces, projectMatrix, weightedTrain] = performStablePCA(trainFaces, meanFace, pcaRatio)
% 稳定的PCA降维，确保特征向量方向一致

% 计算协方差矩阵
C = trainFaces' * trainFaces;
[V, D] = eig(C);

% 按特征值排序
eigenValues = diag(D);
[eigenValues, sortIdx] = sort(eigenValues, 'descend');
V = V(:, sortIdx);

% 选择主成分
cumSum = cumsum(eigenValues) / sum(eigenValues);
numComponents = find(cumSum >= pcaRatio, 1);
if isempty(numComponents)
    numComponents = length(eigenValues);
end

selectedV = V(:, 1:numComponents);

% 计算特征脸（确保方向一致性）
eigenFaces = trainFaces * selectedV;

% 标准化特征脸并确保符号一致
for i = 1:size(eigenFaces, 2)
    eigenFace = eigenFaces(:, i);
    % 确保特征脸的第一个非零元素为正，保证符号一致性
    firstNonZero = find(abs(eigenFace) > 1e-10, 1);
    if ~isempty(firstNonZero) && eigenFace(firstNonZero) < 0
        eigenFaces(:, i) = -eigenFace;
        selectedV(:, i) = -selectedV(:, i);
    end
    
    % 标准化
    norm_val = norm(eigenFaces(:, i));
    if norm_val > 0
        eigenFaces(:, i) = eigenFaces(:, i) / norm_val;
    end
end

% 投影矩阵
projectMatrix = eigenFaces;

% 投影训练数据
weightedTrain = projectMatrix' * trainFaces;

end

function [svmModel, bayesianResults] = trainBayesianSVM(X_train, Y_train, bayesianConfig)
% 使用贝叶斯优化训练SVM模型

fprintf('    正在进行贝叶斯超参数优化...\n');
fprintf('    最大评估次数: %d\n', bayesianConfig.maxEvaluations);
fprintf('    时间限制: %d秒\n', bayesianConfig.maxTime);

% 设置交叉验证分区
cv = cvpartition(Y_train, 'KFold', bayesianConfig.kfold);

% 配置超参数优化选项
optimizationOptions = struct(...
    'AcquisitionFunctionName', bayesianConfig.acquisitionFunction, ...
    'MaxObjectiveEvaluations', bayesianConfig.maxEvaluations, ...
    'MaxTime', bayesianConfig.maxTime, ...
    'CVPartition', cv, ...
    'UseParallel', bayesianConfig.useParallel, ...
    'Verbose', 1, ...
    'ShowPlots', false, ...
    'Repartition', false, ...  % 使用固定的CV分区
    'OutputFcn', @outputFunction);

% 检查MATLAB版本兼容性
if verLessThan('matlab', '9.7')  % R2019b之前的版本
    fprintf('    注意: 检测到较早的MATLAB版本，使用兼容模式\n');
    optimizationOptions = rmfield(optimizationOptions, 'OutputFcn');
end

try
    % 执行贝叶斯优化
    [svmModel, fitInfo, hyperparameterOptimizationResults] = fitcsvm(X_train, Y_train, ...
        'KernelFunction', 'rbf', ...
        'Standardize', true, ...
        'OptimizeHyperparameters', {'KernelScale', 'BoxConstraint'}, ...
        'HyperparameterOptimizationOptions', optimizationOptions);
    
    % 提取优化结果
    bayesianResults = struct();
    bayesianResults.bestParams = struct();
    bayesianResults.bestParams.KernelScale = svmModel.KernelParameters.Scale;
    bayesianResults.bestParams.BoxConstraint = svmModel.BoxConstraints(1);
    bayesianResults.bestObjective = hyperparameterOptimizationResults.MinObjective;
    bayesianResults.numEvaluations = hyperparameterOptimizationResults.NumObjectiveEvaluations;
    bayesianResults.optimizationResults = hyperparameterOptimizationResults;
    
    fprintf('    ✓ 贝叶斯优化成功完成\n');
    fprintf('    最优KernelScale: %.6f\n', bayesianResults.bestParams.KernelScale);
    fprintf('    最优BoxConstraint: %.6f\n', bayesianResults.bestParams.BoxConstraint);
    
catch ME
    fprintf('    ⚠ 贝叶斯优化遇到问题，使用备用网格搜索方法\n');
    fprintf('    错误信息: %s\n', ME.message);
    
    % 备用方案：改进的网格搜索
    [svmModel, bayesianResults] = fallbackGridSearch(X_train, Y_train, cv);
end

end

function [svmModel, results] = fallbackGridSearch(X_train, Y_train, cv)
% 备用的改进网格搜索方法

fprintf('    执行备用网格搜索优化...\n');

% 更精细的网格搜索
kernelScales = logspace(-2, 1, 10);  % 0.01 到 10
boxConstraints = logspace(-2, 2, 10); % 0.01 到 100

bestAccuracy = 0;
bestParams = struct('KernelScale', 1.0, 'BoxConstraint', 1.0);
totalCombinations = length(kernelScales) * length(boxConstraints);
currentCombination = 0;

for i = 1:length(kernelScales)
    for j = 1:length(boxConstraints)
        currentCombination = currentCombination + 1;
        
        ks = kernelScales(i);
        bc = boxConstraints(j);
        
        % 交叉验证评估
        accuracies = [];
        for fold = 1:cv.NumTestSets
            trainIdx = training(cv, fold);
            testIdx = test(cv, fold);
            
            tempModel = fitcsvm(X_train(trainIdx, :), Y_train(trainIdx), ...
                'KernelFunction', 'rbf', ...
                'KernelScale', ks, ...
                'BoxConstraint', bc, ...
                'Standardize', true);
            
            predictions = predict(tempModel, X_train(testIdx, :));
            accuracy = mean(predictions == Y_train(testIdx)) * 100;
            accuracies(end+1) = accuracy;
        end
        
        avgAccuracy = mean(accuracies);
        if avgAccuracy > bestAccuracy
            bestAccuracy = avgAccuracy;
            bestParams.KernelScale = ks;
            bestParams.BoxConstraint = bc;
        end
        
        % 显示进度
        if mod(currentCombination, 20) == 0 || currentCombination == totalCombinations
            fprintf('    进度: %d/%d (%.1f%%)\n', currentCombination, totalCombinations, ...
                    100*currentCombination/totalCombinations);
        end
    end
end

% 使用最优参数训练最终模型
svmModel = fitcsvm(X_train, Y_train, ...
    'KernelFunction', 'rbf', ...
    'KernelScale', bestParams.KernelScale, ...
    'BoxConstraint', bestParams.BoxConstraint, ...
    'Standardize', true);

results = struct();
results.bestParams = bestParams;
results.bestObjective = 1 - bestAccuracy/100;  % 转换为损失
results.numEvaluations = totalCombinations;

fprintf('    ✓ 网格搜索完成\n');
fprintf('    最优KernelScale: %.6f\n', bestParams.KernelScale);
fprintf('    最优BoxConstraint: %.6f\n', bestParams.BoxConstraint);
fprintf('    最优准确率: %.2f%%\n', bestAccuracy);

end

function stop = outputFunction(results, state)
% 贝叶斯优化过程的输出函数

stop = false;
if strcmp(state, 'iteration')
    fprintf('    评估 %d: 目标值 = %.6f\n', results.ObjectiveEvaluations, results.MinObjective);
end

end

function accuracy = testSVMAccuracy(svmModel, X_test, Y_test)
% 测试SVM模型准确率

predictions = predict(svmModel, X_test);
accuracy = mean(predictions == Y_test) * 100;

end

function cvAccuracy = performIndependentCV(X, Y, bestParams)
% 使用最优参数执行独立的交叉验证

cv = cvpartition(Y, 'KFold', 5);
accuracies = [];

for fold = 1:5
    trainIdx = training(cv, fold);
    testIdx = test(cv, fold);
    
    model = fitcsvm(X(trainIdx, :), Y(trainIdx), ...
        'KernelFunction', 'rbf', ...
        'KernelScale', bestParams.KernelScale, ...
        'BoxConstraint', bestParams.BoxConstraint, ...
        'Standardize', true);
    
    predictions = predict(model, X(testIdx, :));
    accuracy = mean(predictions == Y(testIdx)) * 100;
    accuracies(end+1) = accuracy;
end

cvAccuracy = mean(accuracies);

end 